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We investigate the quantum many-body instabilities of the extended Hubbard model for spinless 
fermions on the honeycomb lattice with repulsive nearest-neighbor and 2nd nearest-neighbor density- 
density interactions. Recent exact diagonalization and infinite density matrix renormalization group 
results suggest that a putative topological Mott insulator phase driven by the 2nd nearest-neighbor 
repulsion is suppressed, while other numerically exact approaches support the topological Mott 
insulator scenario. In the present work, we employ the functional renormalization group (fRG) for 
correlated fermionic systems. Our fRG results hint at a strong suppression of the scattering processes 
stabilizing the topological Mott insulator. From analyzing the effects of fermionic fluctuations, we 
obtain a phase diagram which is the result of the competition of various charge ordering instabilities. 

PACS numbers: 


I. INTRODUCTION 

The extended Hubbard model for spinless fermions 
at half-filling is possibly the simplest itinerant fermion 
model with interactions to be put on the honeycomb lat¬ 
tice. Yet, it features interesting interaction driven effects 
as, e.g., the possible realization of an interaction induced 
topological Mott insulatoi^ with a symmetry-protected 
chiral edge state. 

The effect of a repulsive nearest-neighbor interaction 
Vi is comparatively well-understood: Beyond a criti¬ 
cal coupling strength Vi destabilizes the semi-metallic 
regime and the system undergoes a direct and continu¬ 
ous strong coupling quantum phase transition (QPT) to 
a fully gapped charge-density wave state (CDW). The 
critical exponents^^ obtained for this quantum critical 
point (QCP) show that it falls into the Gross-Nevei]^ 
universality class’*'. The QPT corresponding to CDW or¬ 
der was investigated by a variety of methods, ranging 
from a Dyson-Schwinger equation based analysis^ to 
non-perturbative renormalization group calculationJ^HH 
and recent Quantum Monte CarlcP^and newl y dev eloped 
Majorana Quantum Monte Carlo simulationsHUll, 

Considering only a repulsive 2nd nearest-neighbor in¬ 
teraction V 2 , the situation presently ap pear s to be far 
from clear. Previous mean-field theories^EMUl found a 
direct and continuous, strong coupling QPT into a sym¬ 
metry broken phase with a complex bond-order (BO) pa¬ 
rameter Xij^ where ij denotes a 2nd nearest-neighbor 


^ We note that the critical exponents of Gross-Neveu type quan¬ 
tum field theories depend in fact on flavor number and also the 
number of spinor components. Here, we are dealing with two fla¬ 
vors of 4-component Dirac fermions, corresponding to a reducible 
representation of the Dirac algebra in 2+1 spacetime dimensions. 
This universality class is also referred to as Gross-Neveu-Yukawa 
with ^2 order parameter or chiral Ising universality clasi^l^^. 


pair. On a mean-field level, a finite dimerization Xij 
on a 2nd nearest-neighbor bond renormalizes the bare 
hopping matrix and breaks inversion symmetry V (real 
part of Xij) and/or time-reversal symmetry T (imagi¬ 
nary part of Xij)' was founcff^ that the mean-field 
ground-state energy is minimized by the solution break¬ 
ing time-reversal symmetry only. This peculiar state 
would be an interaction induced realization of the Hal¬ 
dane modeP^. The latter is a topologically non-trivial 
model for non-interacting spinless fermions on the hon¬ 
eycomb with nearest and 2nd nearest-neighbor hopping. 
The breaking of time-reversal symmetry is due to a par¬ 
ticular flux configuration, which does not, however, corre¬ 
spond to a finite homogeneous magnetic field penetrating 
the system. The two band model features Chern numbers 
±1 and, correspondingly, for a full bulk gap a quantized 
Hall conductivity, which is carried by a gapless chiral 
edge state at the sample boundary. This state is therefore 
referred to as a quantum anomalous Hall state: it is char¬ 
acterized by a topologically protected chiral edge state, 
but as opposed to the integer quantum Hall state, does 
not require a net magnetic flux through the sample. In 
the sense of the Alt land-Zirnbauer classification of topo¬ 
logical insulators and topological superconductors, the 
Haldane model falls into the unitary symmetry class A^^. 
Therefore, in Ref. [U the mean-field state for the inter¬ 
acting spinless fermion model on the honeycomb lattice 
was dubbed topological Mott insulator. 

Numerical studies of spinless fermions on the honey¬ 
comb and the spinless 7r-flux model on the square lattice 
report the absence of an interaction induced quantum 
anomalous Hall (QAH) phase. In Refs. l2Tll24l the obser¬ 
vation of a direct transition from the semi-metallic (SM) 
phase to a modulated charge density wave phase driven 
by repulsive 2nd nearest neighbor density-density inter¬ 
action was reported. In this work, we follow Ref. I25I26I 
and call this particular type of modulated charge den¬ 
sity wave phase a three-sublattice charge density wave, 
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CDWs for short*. In Ref. [521 authors infer from 
cluster perturbation theory that while QAH correlations 
exist for small clusters, the QAH state ceases to be the 
ground-state for increasing cluster size. The lack of a 
long-range ordered QAH state is further attributed to 
the vanishing density of states at half-filling, which leads 
to an insufficient energy-gain from the formation of a 
QAH dimerization pattern. Another exact diagonaliza- 
tion stud}^^ finds an intermittent Kekule dimerization 
phase, sandwiched between SM and modulated charge 
density phase, but with the QAH state also completely 
squeezed out of the phase diagram. The Kekule dimeriza¬ 
tion phase was recently corroborated beyond mean-field 
theory in the strong-coupling regime by ED studies on 
large cluster^^ and infinite density matrix renormaliza¬ 
tion group (iDMRG)P^. These works further revealed 
charge-ordered ground-states not seen in previous mean- 
field phase diagrams. 

These results rep resen t drastic revisions of the mean- 
field phase diagram^^IMIl] for spinless fermions at half¬ 
filling with Dirac-type low-energy excitations. These 
findings might further indicate the presence of strong 
fermionic and possibly collective bosonic quantum fluc¬ 
tuations close to the would-be transition from the semi- 
metallic state to the QAH state. The precise nature 
of these strong quantum fluctuations, however, remains 
somewhat elusive. The SM to QAH mean-field transi¬ 
tion does not break a continuous symmetr}^, ruling out 
the reduction of possible ordering tendencies in a given 
channel due to the backaction of gapless collective de¬ 
grees of freedom. But the ground-state reported by ex¬ 
act diagonalization for Vi = 0 , V 2 > 0 comes with a huge 
degeneracy in the classical limit. On the quantum level, 
there is, however, only a small finite subset of degenerate 
ground states with different charge arrangements. Since 
only discrete symmetries are broken, also in the case of 
the CDW 3 , we do not expect a strong influence of collec¬ 
tive fluctuations. 

We note, however, that other numerically exact and 
variational approaches employed in Ref. [28] do find sup¬ 
port for the interaction induced QAH state. The differ¬ 
ences between Refs. I21|22l and [28] conce rn system sizes 
and boundary conditions: periodicP^^^ vs. operP^. In 
our functional renormalization group approach to the 
problem, we work with an infinite system with periodic 
boundary conditions. 

In the following, we will present our results obtained 
within the functional renormalization group (fRG) frame¬ 
work for correlated fermions. One of our main results 
in the present work is that the unbiased inclusion of 
fermionic fluctuations with a refined resolution of mo¬ 
mentum space is sufficient, to obtain the GDW 3 with 
a finite wavevector transfer Q as the leading instability 


* Note that in Refs. fT8l2H424l27l the abbreviation for this charge- 
modulated insulating state is CM(s). 



FIG. 1: Left panel: Energy bands for spinless fermions on the 
bipartite honeycomb lattice with nearest-neighbor hopping 
amplitude t. For the system at half-filling, the valence band 
(red) is completely filled while the conduction band (blue) 
is completely empty. The dispersion is approximately linear 
around the K and K' points. Right panel: Lattice geometry 
and bond dimerization Xij ? cf. Ref. [T] White disks corre¬ 
spond to sites of the A sublattice, black disks correspond to 
those of the B sublattice. Red arrows correspond to nearest- 
neighbor vectors (5i, fc, S 3 , and blue arrows correspond to 2nd 
nearest-neighbor vectors Ai, A 2 , A 3 . Thin black lines (solid: 
A sublattice, dashed: B sublattice) visualize the bond order 
Xij • 


from our renormalization group flows, as the interaction 
strength V 2 is increased to destabilize the semi-metal. We 
thus provide further evidence that the sought after topo¬ 
logical Mott insulator phase for spinless fermions on the 
honeycomb is destroyed by competing fermionic fluctu¬ 
ations in the particle-hole channel, and is replaced by a 
fully gapped, charge-ordered state with an enlarged unit 
cell. We further investigate the phase diagram for both 
El, E 2 >0. 

The present paper is structured as follows. In Sect.[TI| 
we introduce the Hamiltonian for spinless fermions on 
the honeycomb lattice with repulsive nearest and 2 nd 
nearest-neighbor interactions. In Sect. iml we recap the 
essentials of the fRG method we employ to analyze the 
phase diagram of the model. In Sect. EYl we present 
our results for the phase diagram, as the strengths of 
both nearest and 2 nd nearest-neighbor interactions are 
varied. Sect.|V|is devoted to a detailed investigation of 
the relevant scattering processes which lead to a suppres¬ 
sion of the QAH instability and ultimately favor ordering 
tendencies corresponding to a charge modulated ground 
state (CDW 3 ). 


II. THE MODEL 


The Hamiltonian H of the spinless fermion model on 
the honeycomb is decomposed into a hopping part, Hq, 
and an interaction part iLint, 


F^ = F^0+^int. 


( 1 ) 
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In the following, we only consider nearest-neighbor hop¬ 
ping with a real hopping amplitude t, 

Ho = (c|cj + cjci) +nY2 

{ij) * 

The interaction term reads 

i^int — hi ^ ^ TliTlj T V 2 ^ ^ TliTlj. ( 3 ) 

(hj) ((id» 

Here, the operator c\ creates a spinless fermion at lattice 
site i. The symbol denotes a sum over nearest- 

neighbor pairs of lattice sites connected by the vectors 5i , 
i = 1,2,3, where each pair is counted only once. Analo¬ 
gously, j)) denotes a sum over 2nd nearest-neighbor 

pairs connected by the vectors A^, i = 1 ,...,6 , where 
again each pair is counted only once. See Fig. for a 
depiction of the lattice geometry. 

Due to the two sublattices with inequivalent nearest- 
neighbor sites, the hopping Hamiltonian gives rise to a 
two-band model. The half-filling condition, which we will 
enforce througout the rest of the paper, simply trans¬ 
lates to a vanishing chemical potential, /i = 0. For van¬ 
ishing interactions, the simple nearest-neighbor hopping 
leads to a band structure analogous to that of the tight- 
binding model for graphene. The Fermi energy is poised 
at the Dirac points iF, K'. The fermionic low-energy 
excitations feature a linear dispersion around these high- 
symmetry points in the Brillouin zone, cf. Fig. where 
the valence and conduction bands touch. Correspond¬ 
ingly, the non-interacting single-particle density of states 
vanishes, as the energy of the non-interacting fermionic 
degrees of freedom approaches the Fermi energy. With¬ 
out interactions, the model can be understood as a 
semi-metal. The vanishing of the non-interacting single¬ 
particle density of states renders the semi-metallic state 
stable with respect to correlation effects, such as the 
opening of a gap in the single-particle spectral function, 
or possible symmetry breaking quantum phase transi¬ 
tions. 

We note that for spinless fermions, the half-filling con¬ 
dition implies one fermion per 2-atom unit cell. The in¬ 
teraction induced quantum anomalous Hall state char¬ 
acterized by a complex bond-order field Xij = 
described in good detail in Ref. [H 

III. FUNCTIONAL RENORMALIZATION 
GROUP ESSENTIALS 

In this work, we employ a functional renormalization 
group approach for the one-particle-irreducible (IPI) ver- 
ticeP^ with an energy cutoff. For recent reviews on the 
fRG method, see Refs. [30]-[32| The fRG calculation is 
most economically performed in the band-basis, which 
diagonalizes the quadratic part Hq of the fermion Hamil¬ 
tonian. The fRG equations can be derived from the func¬ 
tional integral representation of the fermionic partition 


xx=xrrxx 



FIG. 2: Diagrammatic representation of the functional renor¬ 
malization group equation for the scale-dependent 4-point 
vertex V^. The grey filled circle denotes while the black 
dot on the left hand side corresponds to the scale-derivative 
Internal lines correspond to the fermion propagator 
, while slashed lines correspond to the single-scale propaga¬ 
tor S^. Diagrams with slashed and unslashed lines exchanged 
are not shown for brevity. The first diagram on the right 
hand side corresponds to fermion fluctuations in the particle- 
particle channel. The remaining two diagrams encode the 
fluctuations in the particle-hole channel. 


function. The corresponding Matsubara action of the 
model Eq. § can be compactly written as 

+ v['ip,'ip]. (4) 

In the band representation, where the fermion fields cor¬ 
respond to the eigenstates of the hopping Hamiltonian, 
the bilinear part becomes 

{xp, Gq^xP) = y] xpb{k) - eb(fc)) xpb{k), (5) 

b,k 

where k = (icjn: k) and the fermion fields are labeled by 
the band index h = v^c {v: valence, c: conduction band). 
The interaction functional U['0,'0] reads 

ybP,'>P] = ^ E Hi,62,63,64(^l>^2,fc3,fc4) X 

'>Pbi{ki)xpb2{k2)xpb3{k3)xpbAk4)- ( 6 ) 

The properly antisymmetrized coupling function 
^ 61 , 62 , 63 , 64 (^ 15 ^ 5 ^ 35 ^ 4 ) obtained from the interac¬ 
tion Hamiltonian Eq. (§ by substituting operators 
by Grassmann fields and going from the tight-binding 
representation to the band-representation. In the 
process, the coupling function picks up an additional 
momentum dependence - so-called orbital make-up (see 
e.g. Ref. I34|) - due to the non-trivial sublattice structure 
of the hopping Hamiltonian. The frequency dependence 
is simple and just encodes the Matsubara frequency 
conservation of the instantaneous interaction. In the IPI 
fRG-scheme, an infrared regulator with energy scale A 
regularizes the bilinear contribution of the bare action in 
the functional integral. The regularized bare propagator 
replaces the bare propagator according to 

Go(^, ico’n, ^ Gq (^, ico’n, ^) = -^ ^^• (7) 

iuJn - ^h\k) 

The cutoff function is designed to suppress the modes 
with absolute value of band energy \eb{k)\ below the scale 
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FIG. 3: The figure illustrates the discretization scheme for 
the momentum dependence for the external legs of the vertex 
function in the first Brillouin zone for a total of N = 144 
patch points (blue points). The momenta ki, /c 2 , fe are pro¬ 
jected onto the patch points, denoted by 7r{k). These are 
chosen to lie on patch rings around the K, K' points. We 
refer to the patch ring closest to K or K' as the first patch 
ring, the next largest radius corresponds to the second patch 
ring, and so on. See Appendix [B] for details of the projection 
scheme. Due to translation invariance, the fourth momentum 
/c 4 is projected onto the patch point dictated by momentum 
conservation. We evaluate the flow equations Eq. for the 
projected vertex function '\/^(7r(^i), 7r(^2), 7r(fe), 7r(/c4)). The 
scale-resolved loop integrations are discretized along the blue 
lines emanating from the K, K' points in radial direction. For 
clarity, we only enumerate the 36 patch points closest to the 
K, K' points. The patch rings further away from the K, K' 
points are labeled analogously. 

A, 

C^[eb{k)] = Qe{\eb{k) \ - A). (8) 

We chose a smoothed step function 0£, where the 
width of the step is characterized by a softening pa¬ 
rameter 5. The modified scale-dependent propagator 
Gq gives rise to a modified bare action. Performing 
the functional integral yields the effective action T^, cf. 
Refs. l29l31l33l which serves as a generating functional of 
scale-dependent IPI vertex functions. 

By acting with d/dA on P^, one obtains an infinite 
hierarchy of coupled flow equations for the IPI vertex 
functions, in close analogy to the infinite tower of cou¬ 
pled Dyson-Schwinger equations. Integrating the flow 
down from some initial scale Aq we smoothly interpolate 
between the bare action and the effective action at en¬ 
ergy A. Correspondingly, we obtain the IPI vertices of 
the effective action at scale A. The resulting flow equa¬ 
tions are valid for systems with U{1) charge symmetry. 
Their recent use enabled insights into correlated electron 
systems without spin-rotational symmetry, such as co r- 
related electron systems in the spin-orbit regim^^^^ or 
in the presence of magnetic order^. 

Here, we limit ourselves to the rather successful stan¬ 
dard truncation of the infinite tower of differential equa¬ 
tions, suited for analyzing instabilities in correlated 


fermion systems. We keep only the flow of the static 
4-point vertex, i.e., the effective instantaneous interac¬ 
tion V^. We neglect both self-energy feedback and the 
flow of higher-order vertices which are generated upon 
running the mode elimination procedure by integrating 
the flow equations. We will further justify this truncation 
in Sect. I VII 

The flow of the 4-point vertex is given by 

= (ppp + 0ph,d ~ 0ph,cr5 (9) 

where the particle-particle bubble 0pp and the crossed 
and direct particle-hole bubbles 0ph,cr and 0ph,d: respec¬ 
tively, are understood as bilinear functionals of the scale- 
dependent vertex function V^. Fig. shows a diagram¬ 
matic representation of the flow equation Eq. <§■ The 
initial condition for the vertex function is given by 
the coupling function ^ 2 , fe, ^ 4 ) entering 

the Matsubara action. For explicit expressions of the 
flow equation Eq. 0, see Appendix g 

The renormalization of the vertex function obtained 
by solving Eq. © numerically within a so-called patch¬ 
ing scheme yields valuable information about the low- 
energy properties of the model in question and indicates 
instabilities of metallic or semi-metallic phases towards 
symmetry-broken ground-states by a flow to strong cou¬ 
pling. The details of the patching discretization can be 
found in Fig. and Appendix |B| 

From the evolving pronounced momentum structure 
one can then infer the leading ordering tendencies. Since 
we find that the occurring instabilities are most easily 
interpreted in terms of the sublattice basis, after running 
the fRG flow, we finally transform the vertex function 
from the band to the sublattice basis, with the sublattice 
index o = A^B. 

IV. INSTABILITY ANALYSIS 

We solve the fRG flow Eq. (|^ for the effective inter¬ 
action numerically, within the approximations de¬ 
scribed in Sect. in We sweep through values of the cou¬ 
plings Vi and V 2 in the parameter range Vi/t, V 2 /t G 
[0,4]. We note, however, that in the region in parameter 
space where either Vi > 2t or V 2 ^ 2t, the largest com¬ 
ponent of the vertex function fulfills max > 2D with 
the bandwidth D = 6t. In this region, our weak coupling 
truncation of the fRG hierarchy is rendered unreliable. 

In the case of a flow to strong coupling, only certain 
components of the momentum dependent vertex will di¬ 
verge. The divergent momentum pattern, along with 
the corresponding sublattice configurations of the vertex 
function, yield an effective interaction functional. This in 
turn can be translated into an effective low-energy Hamil¬ 
tonian, due to the static approximation for the vertex 
function. From the mean-field ground-state of the effec¬ 
tive Hamiltonian, we obtain a tentative phase diagram 
for the spinless fermion model on the honeycomb in the 
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FIG. 4: Phase diagram for spinless fermions on the honey¬ 
comb lattice at half-filling in the plane of repulsive nearest- 
neighbor and 2 nd nearest-neighbor density-density interac¬ 
tions Vi and V 2 , respectively. The results were obtained with 
a iV = 144 patching-scheme, see Fig. The leading insta¬ 
bility for Vi = 0, V 2 > is the three-sublattice charge 
density wave (CDW 3 ). For Vi > Vi^c, V 2 = 0, the leading 
instability is the conventional charge density wave (CDW) 
with charge imbalance between the A and B sublattices. The 
semi-metallic (SM) state shows an extended region of stabil¬ 
ity. The color-scale corresponds to the critical scales Ac in 
units of the nearest-neighbor hopping t. For very large in¬ 
teractions, the quantum anomalous Hall (QAH) instability 
emerges in the phase diagram. This might be an artefact due 
to the breakdown of the weak coupling approximation. In the 
region marked with N/A, where the three different charge or¬ 
dering instabilities meet, we observe strong K-K scattering. 
A clear identification of the leading instability is this region 
was, however, not possible. 



FIG. 5: Upper panel: Vertex structure in sublattice represen¬ 
tation for divergent GDW correlations for Ui = 1.2t, U 2 = 0 
at the critical scale Ac ~ 0.26t computed in V = 144 patch¬ 
ing scheme as depicted in Fig. I From left to right: 

The patch number of the patch momen- 






ABAB^^^ABBA' 

turn 7 r(A;i) is on the vertical, the patch number of 7 r(A: 2 ) on 
the horizontal axis, while TT(k^) is fixed to the second patch 
(see. Fig. 1^. For clarity, we only display the vertex for the 
36 patches closest to the K, K' points. The divergent mo¬ 
mentum structure can be translated to an effective Hamilto¬ 
nian Eq. (10) favoring a GDW ground state. Lower panel: 
Vertex structure in sublattice representation for divergent 
GDW 3 correlations for Vi = 0 , U 2 = 1 . 2 t at the critical scale 
Ac = 0.70t computed in V = 144 patching scheme as depicted 
in Fig.[^ with the same conventions as above. The divergent 
momentum structure corresponds to enhanced scattering be¬ 
tween fermions with a momentum transfer Q ^ K — K'. This 
translates to the effective Hamiltonian Eq. (11) with a GDW 3 
ground state. The faint checkerboard-like pattern in the back¬ 
ground resembles the initial condition of the vertex function. 
During a portion of the flow, this pattern is also enhanced as 
a whole. 


^ 2 ) parameter space. Our results are collected in the 
phase diagram depicted in Fig.[^ along with estimates for 
the typical energy scales, where long-range order starts 
building up. These are inferred from the renormalization 
group scale where the vertex starts to diverge. We denote 
the corresponding critical scale as Ac- The values for Ac 
can be regarded as estimates for critical temperatures. 

Semi-metallic region. We find an extended region of 
stability for the semi-metallic state. In this regime, the 
flow remains regular and no signs of a flow to strong 
coupling appear. Beyond certain values of Vi, U 2 , the 
semi-metallic state shows instabilities in the particle-hole 
channel toward various charge-ordered states. 

CDW instability. Let us consider the case U 2 = 0 first. 
Beyond a critical coupling Vi^cjt ~ 0-6, the SM state is 
destabilized, and we encounter a CDW instability. The 
diverging momentum structure for the patch-points clos¬ 
est to iF, K' is depicted in the upper row of Fig. The 
CDW instability is driven by scattering of states with 
zero momentum transfer. There is no modulation of the 
scattering amplitude, as the momentum values of the in¬ 
coming and outgoing states sweep through the represen¬ 
tative patch momenta. The effective interaction Hamil¬ 


tonian close to the CDW instability extracted from the 
numerical data for the renormalized vertex reads 

^ 0 , 0 ' eoeo' iV2=o (10) 

o,o' 

where Vo,o' > 0, and Af is the number of unit cells. 

The prefactors ca = +1, = —1 capture the sublat¬ 

tice modulation. The Hamiltonian obviously factorizes 
into a product of density operators k 

for g = 0. This translates into an infinitely ranged 
density-density interaction on the lattice. The sub¬ 
lattice modulation prefers the occupancy of either the 
A or the B sublattice through both attractive intra¬ 
sublattice components ^ ^ 

as well as a repulsive inter-sublattice component ^ 

N^ NN. Since the effective low- 

energy Hamiltonian (unintegrated part of valence and 
conduction bands + effective interaction) conserves par¬ 
ticle number, the resulting ground-state will have either 
the A OT B sublattice occupied/empty. This corresponds 
to the spontaneous breaking of the discrete sublattice 
symmetry, see Fig. for the resulting sublattice occu- 
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pat ion. We note that the diverging momentum pattern 
with momenta on the 1st patch rings closest to the K, K' 
points appears slightly blurred, because in our discretiza¬ 
tion the K-K scattering is represented by all momenta on 
the first patch rings. The momentum pattern on higher 
patch rings, however, becomes very sharp. 

CDWs instability. We now consider the case Vi = 0. 
For strong 2nd neighbor coupling, the semi-metallic state 
turns unstable at a critical coupling of 1 ^ 2 ,c ~ IM. It does 
not, however, correspond to a QAH instability, as o ne 
would expect based on previous mean-field result 
The emerging instability has a rather peculiar momen¬ 
tum structure. We show the diverging momentum struc¬ 
ture in the lower row of Fig. Also here, the momen¬ 
tum structure for momenta on patch-rings further away 
from the K, K' points is very sharp. We find a striking 
Q 7 ^ 0 signature, indicating the tendency toward the for¬ 
mation of unit-cell enlarging order. The intra-sublattice 
component of the interaction vertex, for exam¬ 
ple, is dominated by outgoing momenta /ci, /c 2 = + Q 

(vertical feature) and = ^3 + Q, ^2 (horizontal fea¬ 
ture). This, and also the intermittent amplitude of the 
interaction vertex along these features, correspond to en¬ 
hanced K-K' scattering with a momentum transfer of 
Q K — K'. With the same definitions as above, we 
extract the effective interaction Hamiltonian as 


^eff ~ o'^0^0' 




-Q Q 


«)■(“) 


Upon transforming to real space, one arrives again at an 
infinitely ranged interaction. But in the present case, 
the amplitude shows a modulation with the wavevector 
Q. An analogous effective interaction was already ob¬ 
tained for an e xtended Hubbard model on honeycomb 
bi- and trilayerP^J^^^^. Since the interaction in Eq. 0 
factorizes into a product of fermion bilinears and is 
infinitely ranged, mean-field theory is expected to yield 
reasonably accurate information about the ground-state. 
It turns out that in self-consistent mean-field approach, 
the energy is minimized by a finite complex order param¬ 
eter (A^q) = e^AoC^^, described by its amplitude and 
phase a. As observed already in Ref. [25l this gives rise 
to a density modulation ^ cos (q-r + on the hon¬ 
eycomb lattice, and a concomitant 6 -atom unit cell, cf. 
Fig.§ This implies that each sublattice A, B is broken 
up into three new sublattices Ai, A 2 , A 3 and Hi, H 2 , 
H 3 . For this reason, the notion of a CDW 3 was coined 
in Ref. |25| for this three-sublattice charge-density wave. 
The phase parameter a describes the redistribution of 
charge in each of the emergent sublattices, while the av¬ 
erage remains constant upon changing a. The ground- 
state energy is minimized for a = nn/S with n integer. 
A state with {N^) ^ 0 and a = n 7 r /3 opens a charge 
gap, in agreement with previous resultj^l ^^ l ^^ i 

We note, however, that not only the components of the 
vertex function with momenta close to the K^ K' points 
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FIG. 6 : Density modulations of CDW and CDW 3 patterns. 
White circles denote lattice sites with decreased density, while 
black circles correspond to sites with increased density. The 
radius is a measure for the local deviation from the average 
charge density. Left: The conventional CDW pattern pre¬ 
ferring the occupation of e.g. the B sublattice. Right: An 
example of a CDW 3 pattern with a = tt/S. Around one 
plaquette of the honeycomb, the charge modulation on each 
site differs, corresponding to the breaking of each sublattice 
into three new sublattices. Looking at the plaquette in the 
lower right corner, in counter clockwise direction, we obtain a 
charge modulation of (—(5, -\-S, —S, -h(5, — A, +A), where S and 
A refer to the local charge imbalance and A > S. 


flow to large values. Also the momenta further away from 
the BZ corners grow during the flow, and show strong 
Q 0 signatures. In fact, we find that the inclusion of 
these momenta is crucial in obtaining the CDW 3 as the 
leading instability close to the semi-metallic regime. 

The QAH instability and its Hamiltonian will be dis¬ 
cussed in Sect. |V| We note that for very large U 2 , we 
observe flows with both strong CDW 3 and QAH signa¬ 
tures. For increasing V 2 , the QAH features will even¬ 
tually dominate, at least on the first and second patch 
ring. Since the fRG in the present truncation is a weak- 
coupling method, the appearance of the QAH instability 
in this parameter regime is possibly an artefact of the 
breakdown of the weak-coupling approximation. 


V. SUPPRESSION OF INTERACTION 
INDUCED TOPOLOGICAL INSULATOR PHASE 

Within previous standard patching schemes for the 
analysis of Fermi surface instabilities^ in correlated 
fermion systems, the momentum dependence of the ver¬ 
tex function was typically projected to patch momenta lo¬ 
cated on the Fermi surface. For nodal fermionic systems, 
where the non-interacting degrees of freedom are charac¬ 
terized by a vanishing density of states at the Fermi level, 
this approximation might not capture all relevant contri¬ 
butions to the renormalization of the interaction vertex. 
In the following, we restrict our attention to Ui = 0, 
U 2 >0. 

For the spinless fermion model on the honeycomb lat¬ 
tice, the simplest patching scheme resolves only the mo¬ 
menta closest to the A, K' points, i.e., only the first 
patch ring is taken into account, cf. Fig. This scheme 
resolves the non-trivial angular dependence of orbital 
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FIG. 7: Vertex structure in sublattice representation for di¬ 
vergent QAH correlations for Vi = 0, V 2 = 0.8t close to the 
critical scale Ac ~ 0.31t computed in a V = 36 patching 
scheme. From left to right: ^abab’ ^abba- The 

sharp horizontal and vertical features display the /-wave am¬ 
plitude modulation, see Eq. (12). The number of the patch 
momentum 7 r(ki) is enumerated on the vertical, the number 
of 7 r(/c 2 ) on the horizontal axis, while ^(ks) is fixed to the 
second patch. The checkerboard pattern in the background, 
which resembles the form of the initial condition of the ver¬ 
tex function, also grows to sizable values. The characteristic 
CDW3 signature with a momentum transfer Q ^ K — K' is 
barely visible as a subleading momentum pattern. 


make-up, but neglects the radial momentum dependence 
of the vertex function away from the iF, K' points. Run¬ 
ning the fRG flow (with e.g. N = 36 patches), we obtain 
a momentum signature as depicted in Fig. The effec¬ 
tive interaction Hamiltonian turns out as 


=-^11 Vo,o'eoeo'Sl^^, V ?= o - 


( 12 ) 


with the operator Sj k+q^o k ^ form fac¬ 

tor = sin(\/3/ca:) —2sin(\/3/Ca,/2) cos(3/c^/2). The vari¬ 
ational ground-state of the Hamiltonian Eq. is given 
by a purely imaginary Moving to real space, 

this corresponds to a finite, purely imaginary dimeriza¬ 
tion amplitude Xij on 2nd neighbor bonds, i.e., the QAH 
state. 

Within this approximation to the vertex, the CDW 3 
can be identified as a subleading momentum pattern for 
certain values of V 2 close to h 2 ,c- To analyze the effect of 
shifting the patch momenta to higher patch rings further 
away from the iF, K' points (cf. Fig. on the approxi¬ 
mation, we keep N = 36 patches, but move the position 
of the representative patch momenta closer to the F point 
and run the fRG flow for the vertex function projected 
onto the new set of patch momenta. The critical scales as 
a function of V 2 corresponding to these (arguably unphys¬ 
ical approximations) are collected in Fig. We observe 
that as we move the representative momenta to higher 
patch rings, the critical scales tend to move down, and 
correspondingly, the value of the critical coupling shifts to 
larger values. Eventually, the CDW 3 takes over the QAH 
instability in the full range of V 2 values, as we project the 
momentum dependence of the vertex on the third patch 
ring (counting from the BZ corners). 

Including all N = 144 patch momenta correpsond- 
ing to our best momentum resolution as depicted in 
Fig-i we obtain the blue curve in Fig. for the crit¬ 
ical scales. The CDW 3 is found as the leading instability 




FIG. 8: Red curves: Evolution of critical scales Ac/t from a 
N = 36 patching as a function of V 2 /t for lA = 0, as the 
patch-momenta are shifted closer to the F point. The curves 
are labeled according to the patch ring (cf. Fig. entering 
the approximation. The dashing indicates, when the QAH 
takes over the GDW 3 as the leading instability. While the 
patching with only the second patch ring is susceptible to de¬ 
tecting the GDW 3 in a small parameter range, the patching 
with the third patch ring only yields the GDW 3 over the full 
parameter range. The fourth patch ring yields higher crit¬ 
ical scales than the third again, indicating the importance 
of momenta further away from the BZ corners. Blue curve: 
Gritical scales Ac/t as a function of lA/t for Vi = 0 within 
the N = 144 patching scheme with all momentum configu¬ 
rations of the vertex coupled to each other. Apparently, the 
blue curve is a compromise of the critical scales obtained with 
the one-patch-ring approximation. 


for V 2 < 3.4t. It is justified to say that by including 
momentum configurations with external legs on any of 
the patch rings, and allowing all configurations to talk 
to each other, the QAH instability is in fact suppressed 
over a large portion of parameter space. This is further 
evidenced by the fact that upon inclusion of all patch 
rings, the critical scales actually move down, compared 
to the approximations involving only a single patch ring. 
If no suppression of the QAH signature were at work, it 
would simply diverge, even before the CDW 3 signature 
has cascaded down from the highest to the lowest patch 
ring. 

We further checked our results within another dis¬ 
cretization scheme, which resolves also momenta close 
to the F point. We find the same qualitative behaviour, 
and only small quantitative corrections in the value of 
the critical coupling to larger and the critical scales to 
smaller values. 

We note that the CDW instability remains largely un¬ 
affected upon changing the resolution of our discretiza¬ 
tion. 


VI. CONCLUSIONS & DISCUSSION 

We have to investigated the phase diagram of spinless 
fermions on the honeycomb lattice at half-filling with re¬ 
pulsive nearest and 2 nd nearest-neighbor interactions by 
the functional renormalization group. The fRG repre¬ 
sents a modern implementation of the Wilsonian renor¬ 
malization group idea of successively eliminating high- 
energy degrees of freedom, while tracking the evolution of 
effective couplings in the effective low-energy theory. We 



















compute the flow of the momentum-dependent effective 
interaction (4-point vertex) without self-energy feedback. 
This corresponds to an unbiased resummation of 1-loop 
diagrams in both particle-particle and particle-hole chan¬ 
nels contributing to the effective interaction. Performing 
an instability analysis, we find CDW and CDW 3 insta¬ 
bilities and a suppression of the QAH topological Mott 
insulator. 

Our fRG approach is a suitable tool for the purpose 
of checking the ground-state manifold of the model for 
the presence/absence of the QAH phase in the parame¬ 
ter regime where maxH^° < 2D holds for the following 
reas ons, (i) While mean-field theories with extended unit 
ceii^mSH detect both QAH and modulated CDW phases, 
the absence of QAH in other mea n-field ca lculations and 
in particular exact diagonalizatioii^ ^^ * ^^ * and iDMRGp^ 
points to a competing instability scenario. The fRG, and 
in particular the ‘working-horse’ truncation described in 
Sect. |ml was shown to capture the competition of dif¬ 
ferent ordering tendencies for low-dimensional fermion 
modelJ^ quite reliably, (ii) Dropping self-energy feed¬ 
back and the flow of the 6 -point vertex and higher order 
vertex functions seems justified for the purpose we have 
in mind, as long as the coupling V 2 remains sufficiently 
small. While a partial inclusion of higher-order flow- 
diagrams in the flow of the effective interaction captures 
the effect of collective fluctuations, this kind of truncation 
is only necessary, when one is after realistic gap sizes and 
critical temperatures or an accurate description of critical 
behavior in proximity to a QGP with gapless collective 
excitations. As mentioned already above, in the present 
case, the instability in question (QAH) corresponds to 
breaking of a discrete symmetry, i.e., collective fluctu¬ 
ation effects are expected to play a minor role for the 
competition of QAH and modulated CDW instabilities. 
Thus, the truncation described in Sect. m should suffice 
to correctly identify the fermionic fluctuations driving the 
system critical. 

In summary, we presented strong evidence for the pres¬ 
ence of a direct transition from the SM to the CDW 3 , 
without an intervening topological Mott insulator state. 
For very large V 2 , however, from our fRG results we can¬ 
not rule out that the QAH phase may be stable. Since 
several numerical approacheP^^^ yield phase diagrams 
without any sign of a QAH phase, the dominance of the 
QAH instability over the CDW 3 in our fRG flows for 
very large couplings is most probably caused by the fail¬ 
ure of our truncation to capture strong-coupling effects. 
As the bare interaction strengths grow larger, the rele- 
vant physics can be captured by an Ising-type modeP^^, 
where the hopping acts as a perturbation. At stronger 
coupling, a real-space formulation of the fRG, similar 
in spirit to spin- fRCpl, could provide a more reliable 
starting point in obtaining an accurate phase diagram 
in the strong coupling regime, than our present formula¬ 
tion based on itinerant degrees of freedom. Further, our 
weak coupling study did not detect any signs of a Kekule 
bond-order instability or charge-order instabilities with 


quadrupled unit cells as reported in Refs. I23I24I in the 
parameter range considered. Within the weak-coupling 
regime of the model, we find qualitative agreement with 
recent numerical studiePH^. 

Some other interesting questions also cannot be easily 
accessed within our truncation, such as the order of the 
transition line between GDW and GDW 3 , or the proper¬ 
ties of the corresponding multicritical point. While the 
QGP of the CDW order is well investigated (see Sect.|^, 
much less appears to be known about the QGP corre¬ 
sponding to CDW 3 order and the corresponding effective 
field theory. 

The case away from half-filling was studied e.g. in Ref¬ 
erence [18] by self-consistent mean-field theory. Our fRG 
results for the spinless fermion model on the honeycomb 
at chemical potential /i 7 ^ 0 will be presented in a forth¬ 
coming paper. 

The authors are grateful to S. E. Seidenbecher and 
I. Boettcher for critical comments on the manuscript. 

Appendix A: Flow Equations 

In shortened notation, the bubble expressions are given 
by 

(/.pp = ^V^oLoV\p, (Al) 

</>ph,d = -\v^oLoV\i,,d, (A2) 

<^ph,cr = (A3) 

where the o symbol denotes the channel specific contrac¬ 
tions between loop Kernel L and vertex functions V^. In 
our approximation, the loop Kernel L = is 

built from the bare scale-dependent propagator Gq and 
the so-called single-scale propagator = —d/dKG^. 
The explicit expressions for particle-particle and direct 
particle-hole bubbles read 

</>pp(Ai, 6,6,C4) = ^ n JdriuL{7]2,m,'n3,Vi) X 

(6, Cl, (??4, 7?1,6, ^4 ), 


9^'ph,d(Cl,C2,C3,C4) = -T 11 d7]^L{7]i,r]2,T]3,m) X 

and the crossed particle-hole contribution is given 
through 

0ph,cr(Cl5C2,C3,C4) = 0ph,d(Cl5C2,C4,C3)- (A4) 

Above, we introduced the shorthand Jdr] to represent in¬ 
tegration/summation over loop variables. In the band 
representation, ^ = (6, icj^,/c), 77 = {b 'Since 
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we are interested in ground-state properties, we make a 
static approximation for the vertex function and neglect 
frequency dependence. 

Appendix B: Projection onto Patch Momenta 

The wavevector dependence of the interaction vertex is 
approximated in a so-called A-patch scheme. The Bril- 
louin zone (BZ) is divided into N patches. A given 


wavevector ^ G BZ is projected onto the closest rep¬ 
resentative patch momentum, 7r{k). See Fig. for our 
patching discretization. A single patch is thus composed 
of all the momenta, which have the smallest Euclidean 
distance to the corresponding representative patch mo¬ 
mentum. We then solve Eq. for the projected vertex 
function E^(7r(^i), 7r(^2): 7r(^3), 7r(^4)). We note that 
also depends on band indices b = v^c of outgoing (first 
two arguments) and incoming (last two arguments) legs. 
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